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Abstract 

We demonstrate that geodesies in exact vacuum Kundt gravitational waves may ex- 
hibit a highly complicated behaviour. In fact, as in the previously studied case of non- 
homogeneous pp- waves, for specific choices of the structural function the motion appears 
to be genuinely chaotic. This fact is demonstrated by the fractal method. 
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1 Introduction 



Chaotic motions were already found and described in a number of various spacetimes, in partic- 
ular in the black hole systems, see e.g. [1-13], or in special axisymmetric solutions [14-17]. Some 
time ago we discovered chaotic behaviour of geodesies in specific vacuum pp- waves [18-21]. Al- 
though chaotic interaction of test particles with linearized gravitational waves propagating on 
flat Minkowski background had been studied before [2,6,22-25], this was the first explicit 
demonstration of chaos in exact radiative spacetimes. 

A natural question thus emerged whether a similar complicated behaviour of geodesies 
could be found in other exact solutions of Einstein's equations that represent gravitational 
waves. There exist several important classes of such solutions which are geometrically different 
from pp-wayes, for example the family of Kundt spacetimes with non-expanding wavefronts or 
radiative spacetimes of the Robinson-Trautman type [26] which describe expanding waves. 

In the present work we concentrate on specific vacuum type N spacetimes from the class of 
solutions that admit a non-expanding, shear-free, and twist-free null geodesic congruence which 
was first described in 1961 by Wolfgang Kundt [27-29] (for recent reviews see, e.g., [26,30]). 
These so-called Kundt waves may be interesting from the theoretical point of view because 
all curvature invariants of all orders identically vanish for such spacetimes [31,32]. They may 
thus play an important role in string theory and quantum gravity since there are no quantum 
corrections to all perturbative orders [33]. 

Some aspects of geodesic motions in these spacetimes have already been studied and de- 
scribed. For example, it was shown in [34] that geodesic deviations exhibit a transverse struc- 
ture of relative motions of test particles. In [30] we presented explicit families of timelike, 
null and spacelike geodesies and parallelly propagated frames which enabled us to elucidate 
some global properties of the spacetimes, such as an inherent rotation of the wave-propagation 
direction, or the character of singularities. In particular, we demonstrated that the character- 
istic envelope singularity of the rotated wave-fronts [35] is a (non-scalar) curvature singularity, 
although all scalar invariants of the Riemann tensor vanish there. 

Nevertheless, an open question has remained whether geodesic motion itself exhibits chaotic 
behaviour, analogous to the chaos in pp-waves [18-21]. In the present work we demonstrate 
that this is indeed the case. After a brief review of the Kundt vacuum waves in section [2j we 
show in section [3] by a numerical fractal method that geodesies in these spacetimes are chaotic. 
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2 Kundt vacuum spacetimes of type N 



The metric of the Kundt gravitational waves [26-29] (distinct from the widely known class of 
pp-waves) can be written in the form [30] 

ds 2 = dx 2 + dy 2 - 4x 2 du dv + A(x 2 v 2 + xG) du 2 , (1) 

where the function G(x,y,u) depends on retarded time u and on spatial coordinates x,y 
which span the plane wave-fronts u = const. Due to the presence of a (non-scalar) curvature 
"envelope" singularity at x = 0, we have to restrict ourselves to the region in which x > 0. The 
geodesic equations corresponding to (JXJ) are 

x = -Axvii + (Axv 2 + 2G + 2xG jX )u 2 , (2) 

y = 2xG^ii 2 , (3) 

u = -2-u - 2vix 2 , (4) 
x 



-2-v + Avvu + 2 (— \ xu + 2(—\ yu + 



— \ - 4t> 3 - 4v— u 2 , (5) 

& \ *b / ,x \ / .y V «^ / ,u *^ 

where the dot denotes differentiation with respect to an affine parameter r. In addition, there 
is the normalization condition of four-velocity: 

x 2 + y 2 - Ax 2 iiv + 4(x 2 v 2 + xG) u 2 = e , (6) 

where e = —1, 0, +1 denotes respectively the timelike, null or spacelike character of the geodesic. 

Specific geodesic motion given by equations ([2])-([6]) crucially depends on the particular 
form of the spacetime function G(x,y,u). We concentrate here on type N vacuum Kundt 
spacetimes ([T]). It such a case the function G/0 satisfies the Laplace equation G )XX + G m = 0, 
see [30]. A general solution to this field equation is thus G = Re{g(£,w)}, where g is an 
arbitrary function holomorphic in the argument £ = x + i y. If the function g is linear in £ or 
independent of it, this generates only a flat space. The simplest radiative Kundt spacetimes 
thus arise when g = d(u)(£ — £ s ) n , n = 2, 3, 4, . . . , £ s = const. As in [18, 19] we will consider 
only the waves of a constant profile, i.e. G independent of u, so that d(u) = A = const. We can 
then use the translational symmetry in y to remove the constant y s . We will thus investigate 
geodesies in "polynomial" Kundt spacetimes characterized by the structural function 

G = XRe{(x-x s + iy) n } , (7) 

where A, x s are constants, and n — 2, 3, 4, • • • . Introducing polar coordinates, x — x s = r cos<^, 
y = r sin (p, we may write G = A r n cos(n<^) which is an n-saddle function. 
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The metric ([!]), ^ admits the Killing vector d u whose norm is given by g uu = A(x 2 v 2 + xG). 
These Kundt spacetimes can thus be invariantly divided into regions in which this Killing 
vector is timelike or spacelike, separated by the boundaries on which the vector d u in null. 
This appears when xv 2 is smaller tan, greater than, or equal to —G(x,y), respectively. For a 
fixed value of v this is visualized in figure [T] for n — 6,7, 8: the upper pictures correspond to 
x s = while the lower ones correspond to x s = 1. The Killing vector is timelike in the gray 
regions while it is spacelike in the white regions of the (x, |/)-plane. 




(d) n = 6 (e) n = 7 (f) n = 8 



Figure 1: Regions of different character of the Killing vector d u : in the gray regions it is 
timelike, while in the white regions it is spacelike. On the boundaries the Killing vector is null. 
The upper three pictures correspond to x s = 0, the lower ones correspond to x s — 1. 

3 Demonstration of chaos in polynomial Kundt waves 

The system of non-linear differential equations ([2])-([5]) with Q is quite complicated, so one can 
explicitly find only a few very special analytic solutions, see [30]. In order to obtain the global 
picture of a geodesic motion in the studied spacetimes, it is necessary to perform numerical 
investigations. 
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3.1 The choice of initial data 



In these numerical integrations we make the following choice of initial conditions. In the Kundt 
coordinates (x,y,u,v), the initial position of each geodesic x M (r) is given by 

a^(0) = (x , 2/0, (W 1 ) , (8) 

where 

x = x s + r cos <f , y = r sin (p , (9) 

while the initial velocity is 

— (o) = (-— o — — L + J. + G ( 3 -'o^») N i (10) 

For given x s and r , it is thus a one-parameter family of geodesies parametrized by ip. 

The geometric meaning and justification of such "somewhat peculiar" choice of the initial 
conditions is provided by the transformation from the Kundt coordinates (x, y, u, v) to the 
standard Minkowski coordinates (X, Y, Z, T) for the flat space case when G = 0, namely 

X = x{l + 2uv) , Y = y , 

Z = x[v — u(l + uv)] , T = x[v + u(l + uv)] , (11) 
see [30,35]. We immediately obtain 

X"(0) = (x , ye, 1, 1) , -frr(°) = (°> °> °> X ) ' ( 12 ) 

All these test particles thus start from rest (at the same time T = 1 and plane Z = 1) from 
the position (X, Y) = (x ,y ), i.e. on a drc/e of radius r , centred around X = x s , Y = 0. 
The parameter <£> in ([9]) is the standard polar angle along this circle. In view of the geometric 
interpretation of the Kundt spacetimes presented in [30,35], the initial positions of the geodesic 
test particles are localized on u = 0, which is a half-plane tangent to an expanding cylinder 
X 2 + Z 2 — T 2 = x 2 = 0, Y arbitrary, as illustrated in figure |2j 



The last term in expression (10), which is proportional to G, has been included in order to 
guarantee the consistency of the geodesic equations and the normalization condition (|6j) across 
the wave-surface u = even in a non-flat case. In other words, we assume the continuity 
of the geodesies in the Kundt coordinates (x,y,u,v) on u — 0, and also the continuity of 
the corresponding coordinate velocities, except for an inevitable jump in v equal to (G/x)u. 
The physical interpretation of these assumptions is that, in fact, we investigate geodesic motion 
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(a) G = XRe{(x + iy) 8 } (b) G = ARe{(x — 1 + iy) 8 } 

Figure 2: Initial positions of test particles are localized on a circle of radius r in the planar 
wave-surface u = which is tangent to an expanding cylinder with the vertical axis Y and 
transverse directions X, Z. The particles are at rest in Minkowski space in front of the Kundt 
(shock) wave. The value of the metric function G(x, y) on the initial wave-surface is also 
indicated by various levels of grey (cf figure [T]). 

of free test particles which are initially at rest in a (locally) Minkowski space in front of a shock 
wave. These particles are hit simultaneously by the initial wave-front u — 0, and subsequently 
they start to move under the influence of the Kundt gravitational wave described by Q, @- 

3.2 The geodesies 

Let us now present the numerically obtained results which demonstrate a typical behaviour of 
geodesies in the above "polynomial" Kundt spacetimes. 

In figure [3j we show a family of geodesies in the spacetimes ([TJ, ([7]) with x s = 0, n — 8, 



i.e. G = ARe{(x + iy) 8 }, for initial data (p])-(10), ignoring here the value of the coordinates 



u and v . We have chosen the characteristic radius r = 1 and considered the cases when 
A = 2,4,16,32 which are shown in parts (a), (b), (c), (d). 

We observe that the geodesies either (in a finite proper time) terminate in the envelope 
singularity x = 0, or are unbounded and radially escape to infinity in the (x, ?/)-plane. Interest- 
ingly, they escape through only one of the outcome channels which exactly correspond to the 
axes of the regions in which the Killing vector d u is spacelike, cf the white regions in figures [l] 
and [2] As in the case of non-homogeneous pp-w&ves [18-21], as the test particles escape, their 
frequency of oscillations around the axis of the outcome channel grows to infinity while the 
amplitude of oscillations tends to zero. For small values of the parameter A, the "attractive ef- 
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(a) 



(b) 




0123456 0123456 



(c) (d) 
Figure 3: Trajectories of the test particles for G = A Re{(x + iy) 8 }- 

feet" of the envelope singularity is dominant so that almost all particles approach x = while, 
with a growing value of A, the effect of the gravitational wave becomes stronger and thus more 
and more geodesic particles are forced to escape to infinity through the outcome channels. In 
fact, they are attracted to the curvature singularity at r = oo, see expression ^ which has 
the form G = A r n cos(nip). 

This characteristic behaviour also occurs in the case G = A Re{(x — 1 + i y) 8 } corresponding 
to x s = 1, n = 8, with the radius of the initial circle r = |, illustrated for A = 8, 16, 32, 64 in 
parts (a), (b), (c), (d) of figure |4j respectively. Again, for small values of the parameter A 
all the geodesies end in the envelope singularity x = 0, but with a growing A the effect of 
the polynomial Kundt gravitational wave increases and the outcome channels open. These 
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0.5 1 1.5 2 2.5 3 3.5 4 0.5 1 1.5 2 2.5 3 3.5 4 



(c) (d) 
Figure 4: Trajectories of the test particles for G = A Re{(x — 1 + iy) 8 }. 

coincide with the regions in which the Killing vector d u is spacelike, see figures [T], [2j Although 
the geodesies which move through the channels with decreasing x still terminate at x = (in 
fact, such channels are somewhat "bent" by the presence of the envelope singularity), those 
with a growing x escape to infinity, r — > oo. 

3.3 The fractal structure of geodesic motion 

The main objective of the present paper, however, is to establish the chaotic behaviour of 
geodesies in the Kundt-wave spacetimes. 

Chaos is usually indicated by a highly sensitive dependence of possible outcomes on the 
choice of initial conditions. The standard invariant approach to prove chaotic motion, called 
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a fractal method, was advanced in general relativity, e.g., in the papers [1,3-5]. It starts with 
a definition of several different outcomes, that is "types of ends" of all possible trajectories. 
Subsequently, a large set of initial conditions is evolved numerically until one of the outcome 
states is always reached. Chaos is established if the basin boundaries which separate initial 
conditions leading to different outcomes are fractal. We shall now demonstrate that exactly 
these structures can be observed in the system studied. 

A natural parametrization of the unit circle of initial positions Q, ^ for the geodesies 
(some of which are visualised in figure [4]) is provided by the polar coordinate <p. Obviously, 
for large values of A there are (at most) n outcome channels throughout which the unbounded 
geodesies may approach the curvature singularity at infinite values of r, where n is the power 
of the metric function ([7]). These distinct channels represent possible outcomes of our system, 
and we assign them the symbol j. In the illustrative case n — 8, it takes one of the correspond- 
ing discrete values, j G {—2,-1,0,1,2}. In particular, j = denotes those geodesies which 
approach infinity through the channel centered around the axis y — 0, j = 1 corresponds to 
"upper right" channel, j = 2 corresponds to "upper left" channel, j = — 1 and j = —2 are their 
"lower" counterparts, see part (a) of figure |5j 




X initial data i 



(a) Trajectories in the (x, y)-plane. (b) The function j((p(i)). 

Figure 5: Trajectories of test particles which start from the part of the circle of radius 
r = | centered around (x s ,y s ) = (1,0). The angular variable if changes discretely from 
|| to |: there are 800 particles labeled by i — 1,2, .. . ,800 whose initial positions are given 
by (p(i) = ff + 24uoo^- ^he bold line between i = 350 and i = 450 on the horizontal axis of 
part (b) denotes the region which will be studied in detail in the following figure [6j 

In addition to these five labeled outcome channels, there are three undistinguishable chan- 
nels oriented towards the left. These are influenced by the envelope singularity localized at 
x = so strongly that we must ignore them as outcome channels: we assign them a common 
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superfluous label "—4". To be precise, a geodesic is assigned the outcome value j = —4, once 
the value of its coordinate x drops below 0.5 — this (somewhat artificial but given) boundary 
is denoted by a vertical dashed line in part (a) of figure [5j The dashed circle of radius 1.5 
indicates the value of r at which the labels j = —2, —1,0, 1,2 are assigned to the remaining 
geodesies, i.e. where the outcome channels are discriminated. Note that the function j(ip) 
is the analogue of the scattering function which is usually introduced in the classical chaotic 
scattering problems [36] with unbounded trajectories. 

We observe from the right part (b) of figure [5] that in certain regions the behaviour of 
the function j((f) depends very sensitively on if. As a typical example we present j(ip(i)) 
calculated numerically for 800 geodesies labeled by i whose positions on the initial circle 
r = | (shifted by x s = 1 in the spacetime with A = 512 and n = 8) were parametrized 
by (p G (ff , f ) = (0.4189, 0.5236) such that cp(i) = ff + where i = 1, 2, . . . , 800. We dis- 



play some of the corresponding trajectories in the left part (a) of figure [5] In fact, the basin 
boundaries between the outcomes — indicated by the discontinuities in the function j(ip) - 
appear to be fractal. This is confirmed on the enlarged detail of the part (b) of figure [5j and 
on the enlarged detail of this detail in figures [6| [7] and [8| respectively. 




0.5 1 1.5 



(a) Trajectories in the (x, y)-plane. 



X X Xt 



100 200 300 400 500 600 700 800 

initial data i 

(b) The function j(tp(i)). 



Figure 6: Refinement of the previous figure |HJ The angular variable cp changes from |^|| to 
according to ip(i) = |||| + 18 Q QQQ i, where 2 = 1,2,..., 800. The bold lines on the horizontal axis 
of part (b) depict the two regions studied in detail in figures [7] and [8] (index i running from 
131 to 171, and from 584 to 637, respectively). 



In figure [6] we show such zooming in the interval if G 
we plot 800 geodesies initially located at <p(i 
nification covering the smaller interval (p G 
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(0.4639, 0.4779) where 
In figure [7] there is a deeper mag- 
(0.4662, 0.4669) by 800 geodesies 
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(a) Trajectories in the (x, y)-plane. 



100 200 300 400 500 600 700 800 
initial data i 

(b) The function j(tp(i)). 



Figure 7: Refinement of the previous figure [6j The angular variable <p changes from |§§|§f to 



200637T 
135000 



. There are 800 particles numbered by % which start at ip(i) 



200337T 
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3600000 



I. 




11 



starting on tp{%) = |||||^ + 360 oooo ^' an d similarly in figure [8] we present results for another 800 
geodesies which start in the interval if G (f^jf , 135000 ) = (0-4741, 0.4750) with initial positions 
^(*) = 45000 + 2700000 ^ From these typical pictures we conclude that the basin boundaries sep- 
arating the initial data leading to different discrete output channels j are highly complicated 
and, in fact, appears to have a fractal structure. The dependence of the outcome on the initial 
condition is extremely sensitive. 

It thus seems to be a plausible conjecture that the geodesic motion in the Kundt gravi- 
tational wave spacetimes studied is genuinely chaotic, as in the previously investigated case 
of other spacetimes [1,3-5], or [18-21]. However, to rigorously prove this statement would 
require a demonstration that the basin boundaries are fractal to an arbitrary high level, i.e. 
for differences in the initial data which approach zero. Of course, this cannot be achieved 
solely by numerical integration of geodesies due to the inevitable numerical errors. In the 
case of chaotic motion in non-homogeneous pp-w&ves [18-21] the numerical investigations were 
complemented and supported by analytic methods: in the transverse section spanned by the 
(x, y) spatial coordinates the evolution equations are equivalent to a Hamiltonian system with 
Henon-Heiles-type or more general n-saddle polynomial potentials which are "textbook" ex- 
amples of chaotic systems [37-39]. Unfortunately, in the present case of motion in the Kundt 
spacetimes, no such analytic works are available. The reason is that the system of equations 
(2)-(5) is much more complex. In particular, the motion in the transverse (x, y) section is not 
decoupled from the u and v coordinates. 

We can thus present only the following heuristic argument for the presence of a genuine 
chaos. From our numerical simulations we observe that in the (x, y) section the geodesies are 
asymptotically restricted to the regions in which the Killing vector d u is spacelike, see figure [TJ 
Moreover, the geodesies "bounce" on the corresponding boundaries, as indicated on figure [9] 
which shows the details of figures [5j [6j [7] and [8j Let us note that the boundaries also depend 
on the value of the coordinate v and they therefore evolve with time. In the central region 
these geometric boundaries effectively form walls with (time dependent) concave shapes, and 
it is also natural to expect that the geodesic particles may bounce arbitrarily many times 
between them before they are scattered to one of the outcome channels labeled by j. The 
curved wall segments have a "dispersing effect" on the incident trajectories, and this leads to 
a rapid separation of nearby geodesies after a few bounces. We thus expect that the basin 
boundaries separating the outcomes are fractal to any level (the level corresponding basically 
to the number of bounces performed between the walls in the central region), as in the familiar 
chaotic billiard systems [36]. 
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0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 
(a) Detail of the figure p\ 




0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 
(b) Detail of the figure [6] 




0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 
(c) Detail of the figure [j] 




0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 
(d) Detail of the figure [§] 



Figure 9: Zoom in the (x, y)-plane which shows details of parts (a) of figures pj |6j |7| and [8 
The geodesic test particles effectively bounce on concave boundaries which results in a rapic 
dispersion of the trajectories and their fractal separation. The higher the level of the fractal, 
the greater the number of bounces that a particle undergoes before it escapes to infinity. 
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4 Conclusions 



Some time ago, the first explicit demonstration of chaos in exact radiative spacetimes was given 
for a family of non-homogeneous pp-w&ve solutions which represent the simplest models of exact 
gravitational waves in general relativity [18-21]. In the present paper we have extended these 
results. Using the invariant fractal method we have numerically demonstrated a complicated 
chaotic character of geodesic motion in the class of Kundt spacetimes ([TJ with polynomial 
metric function ([7]). These vacuum spacetimes are of type N, and therefore they also represent 
exact gravitational waves. 

Our work shows that geodesic motion in spacetimes even as simple as pp-w&ves and the 
Kundt waves can be complex for specific but natural choices of the corresponding metric 
functions. This may possibly stimulate future investigation of chaotic behaviour of geodesies 
in other exact radiative spacetimes, such as non-expanding Kundt waves with a cosmological 
constant. 
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